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Abstract 

Molecules like water have vibrational modes with zero point energy well above room temperature. 
As a consequence, classical molecular dynamics simulations of their liquids largely underestimate 
the kinetic energy of the ions, which translates into an underestimation of covalent interatomic dis- 
tances. Zero point effects can be recovered using path integral molecular dynamics simulations, but 
these are computationally expensive, making their combination with ab initio molecular dynamics 
simulations a challenge. As an alternative to path integral methods, from a computationally simple 
perspective, one would envision the design of a thermostat capable of equilibrating and maintaining 
the different vibrational modes at their corresponding zero point temperatures. Recently, Ceriotti 
et al. [Phys. Rev. Lett. 102, 020601 (2009)] introduced a framework to use a custom-tailored 
Langevin equation with correlated-noise that can be used to include quantum fluctuations in classi- 
cal molecular-dynamics simulations. Here we show that it is possible to use the generalized Langevin 
equation with suppressed noise in combination with Nose-Hoover thermostats to achieve an efficient 
zero-point temperature of independent modes. We apply this new thermostat to the molecular dy- 
namics of a flexible water force field. We address the question of whether thermostating each mode 
to its zero point temperature is enough to simulate nuclear quantum effects in water. 

PACS numbers: 
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I. INTRODUCTION 



Understanding how large are zero point nuclear quantum effects (NQE) both in water-"- 
and ice^ is an active area of research. How much the structure of liquid water is dependent 
on the classical treatment of the ionic degrees of freedom in ab initio molecular dynamics 
(MD) simulations is an open question 3,6 . Even if a number of path integral molecular dy- 
namics studies have addressed the issuei"- 13 -^, a definite answer has not yet been provided. 
The problem is subtle, due to the complex nature 7 of the OH-0 hydrogen bond (Hbond) 
in water. It is well known that hydrogen bonded materials show an anti-correlation- be- 
tween the high energy, stretching frequencies and the librational frequencies of the molecules. 
Recently^, we have shown that this anti-correlation is the origin of negative griineisen pa- 
rameters of the high energy vibrational modes in ice. These are large enough to cause an 
anomalous isotope effect in the volume of ice, making the volume per molecule of heavy or 
D2O ice larger than that of normal or H2O ice. This anomaly is not captured by flexible 
and/or polarizable force-fields, due to their underestimation of the anti-correlation effect^. 
Nonetheless, we choose to use in this study the q-TIP4P/F 1 force field. Even if it has been 
shown to fail in the description of the anomalous isotope effect of ice^, it provides a good 
qualitative description of the anharmonicities of all the modes in liquid water. In classical 
MD simulations of force field models, all the modes are equilibrated at a given constant 
temperature. This equipartition of temperature is a classical description of liquid water 
which lacks NQE. Recently, Ceriotti et air-— have shown that the key features of path 
integral molecular dynamics (PIMD) simulations of liquid water can be reproduced using 
custom tailored thermostats based on generalized Langevin dynamics (GLE). In their work, 
they were able to enforce the w-dependent effective temperature T(co) = jjf^ coth 2 j^ T si- 
multaneously on different normal modes, without any explicit knowledge of the vibrational 
spectrum. The tailoring aspect of their thermostat involves complicated optimization to 
independently tune the drift and diffusion parameters of the GLE dynamics. In this work 
we introduce a new thermostating scheme with very few, and easy to tune parameters that 
can equilibrate modes to different temperatures. In our scheme we couple both Nose- Hoover 
(NH) and GLE dynamics to the system. We use GLE kernels that satisfy the fluctuation- 
dissipation (FD) condition which can be derived from a well defined harmonic bath model. 
We suppress the noise term in GLE dynamics by setting the GLE temperature close to 0. 
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In this limit, the dynamics is almost deterministic. The frequency dependent equilibration 
is achieved through the independent tuning of NH and the frequency dependent friction 
profile. Microscopic details of the full dynamics are presented in Sec. (Ill A|) . We sacrifice 
transferability of parameters between different systems in exchange for simplicity in their 
optimization against the known vibrational spectrum of the system. This thermostat acts 
on the system within a deterministic regime and hence our method can be thought as a 
deterministic frequency dependent thermostat or phonostat^. The goals of this study are 
two sided. On the methodological side, after rigorously deriving the thermostat equations, 
we evaluate its performance , by comparing it with PIMD simulations of q-TIP4P/F water. 
In addition, we address the question of competing quantum effect si or competing anhar- 
monicities in water- using a quantified, temperature-dependent approach. To achieve this 
we reformulate the idea of NQE in terms of the zero point energy of individual modes. 

II. NOSE-HOOVER THERMOSTATS IN PRESENCE OF GLE KERNEL 

In this work, we construct a new frequency dependent thermostating scheme. This scheme 
involves use of two thermostats. We use the standard NH chain thermostats which is the 
gold standard of thermostats. To achieve frequency dependent equilibration, we use GLE 
to modify the NH action. Recently, Ceriotti and Parinello*^— developed an extensive 
thermostating scheme based on GLE dynamics. We choose a particular form of GLE dy- 
namics and use it in conjunction with the NH dynamics to enforce frequency dependent 
thermostating. We call this new scheme NGLE (for Nose-GLE) thermostating. 

A. Microscopic derivation of NGLE dynamics 

In this section we start with the microscopic model for NGLE thermostat starting from 
a system-bath coupling model. The full extended Hamiltonian of the system can be written 

as, 

Pi 

Htotai = H sys {—, qi) + H NH (p s , s) + H GLE (pi yXk ,x i)k ) (1) 
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Where, 



3N 2 

E - = E^ + ^-^iv) (2) 

i=l 1 
2 

H NH = ^ + (3N + l)k B T NH log s 

3N 

H GL E = y^#GL£ (3) 

(4) 

s is the parameter that modifies the effective dynamics and enforces constant temperature 
ensemble on the rescaled momenta (-). Q is the NH mass. This rescaled dynamics is also 
coupled to a harmonic bath that enforces generalized Langevin dynamics. This is achieved 
by coupling each system degree of freedom to g harmonic oscillators of mass m and frequency 
Uk- /3 is the coupling strength of the oscillator to the system degree of freedom. Pi tXk is the 
momentum conjugate to the kth oscillator position (index i corresponds to the system 
degree of freedom). Hamilton Jacobi equations for the total dynamics for the system degrees 
of freedom can be written as, 

q * = w? (5) 



Pi 



dV(q) ^ (3 



dqi ^— ' mw 



f3 2 



The dynamics of NH degrees of freedom is given as, 

Ps 

S = Q 

^ V ] (3N + l)k B T NH 

i 1 

We can further simplify NH dynamics in Eq. ((7j) by rescaling time and momenta. We perform 
dt —7- — in Eq. (Q. The resulting equations can be written in terms of the rescaled momenta 



Pn 

" = Q 

N 2 

Pv = [Y.J[-^N+l)k B T NH 



(8) 



Where 77 = log s. 

The dynamics of GLE bath degrees of freedom is given by. 



Pi,x k 



(9) 
(10) 



Pi,x k 



m k 

-f3qi - mulx itk 



The above equations for the GLE bath can be solved exactly and can be substituted in 
Eq. (Q. The resulting system's dynamics in presence of NH chains and the GLE bathM^- 
can be written in the following form, 



where the last term is the NH term that is coupled to the system. The memory kernel K(t) 
has an exact expression in terms of the bath parameters, 



The 'random' force £(t) can also be completely determined in terms of the bath degrees of 
freedom, ((t) is connected to the memory kernel through the FD theorem as {({t)((t')) = 
ksToLE K(t — t') with (((t)) = 0. Note that the temperature enforced by the FD condition 
is different from the NH temperature. We can rewrite Eq. (TT2l) by absorbing the NH term 
into the time integral. 




(11) 



(12) 




(13) 




(14) 



Where we have defined, 



K(t,t') = K(t 



f) + 5{t - £') 



(C(t)C(o)) = Kit) 



(15) 



Q 



Now consider Eq. ( fT4"l) in the case when Tqle — > 




(16) 
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In this case the noise term is completely suppressed and we have a GLE dynamics that 
involves only a nonlocal friction profile and the deterministic NH term that provides constant 
temperature T^h to the system. The conserved energy for this modified dynamics can be 
written as, 

2 

H' = H sys {p, q ) + ^L + (3N + l)k B T NH r) + £ AKE t 

(17) 

The last term in the above equation is the change in kinetic energy for each GLE action 
summed over all past trajectories^ 1 ^. In this work we achieve frequency dependent equili- 
bration as a result of the competition between the NH dynamics and the damping action of 
nonlocal friction profile with suppressed noise. 



B. Delta-like memory kernels 

The friction term in Eq. (|12|) is linear in the system momentum, and the friction coeffi- 
cient K(t) is a simple function of the frequencies uj^ and the coupling constants This 
generalized Langevin equation is exact and its validity is not restricted to small departures 
from thermal equilibrium. Now we consider the case of infinite number of oscillators with 
continuous distribution of frequencies Uk- In this limit the summation in Eq. (1T31) can be 
replaced by an integral with some distribution function (^ — >■ N J dug(oj)). 

K(t) = N J du g(u) (-y=y cos ( w *) ( 18 ) 

Since we have control over the bath degrees of freedom we choose the coupling constant to 
~3m = \Z 7 ^ Aa jv + ~- We are free to choose the distribution of frequencies of the oscillators 
to enforce desired memory kernel on the system. To design delta like memory kernels K(t) 
introduced in Ref. (jl2l ). we choose the frequency distribution function of the oscillators to 



be, 



to 2 



Aoj z + [co — Uq) z 

The above distribution of the continuum oscillator frequencies results in the effective delta 
like friction profile of Ref. (112h which is the essential component of frequency dependent 
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thermostating scheme. The memory kernel obtained is given by, 

Aw 




Alo 2 + (lo — uj y 



+ Au 2 + (u + tu ) 2 ) (20) 
K(t) = ^(^^y-^coscoot (21) 

In the memory kernel defined above, 7 is the friction coefficient (or coupling strength to the 
harmonic bath), ujq is the frequency at which the delta shaped memory kernel has maximum 
friction value or maximum strength of coupling. Au is the width of the friction profile. 
Notice that all the parameters related to the GLE thermostat are completely independent 
of the force field parameters. All we need to know is the position of the peaks of modes 
that we can obtain from the spectral density of the system in consideration. Non-Markovian 
dynamics can be mapped to Markovian dynamics in higher dimensional spaced by adding 
auxiliary momentum degrees of freedom^. The modified higher dimensional dynamics can 
be implemented in the form of following dynamical equations, 

q = 2- (22) 
m 

~fj- A r s j + B ^) ( 23 ) 

Matrices A and B are the drift and diffusion matrices respectively. £ is the uncorrelated 
Markovian noise, and s is the vector of additional momentum degrees of freedom. The drift 
and diffusion matrices may be constrained by the generalized FD theorem 

(A + A T ) = Mk B T GLE BB T (24) 

The matrix A has the form | _ | . One can obtain functional form of the memory 

kernel from the matrix A. 




T 

a ps a 



K{t) = 2a pp 5(t) - a ps e-l*l a aj s (25) 
The drift matrix A for that produces delta like friction profile in Eq. (1201) is given by, 

/ v^W) \f^W)\ 



A 



V-VT(^) --0 Ac J 



(26) 
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C. Simulation Details 



In this section we describe the implementation of NGLE thermostat. We begin by writing 
the classical evolution for the full dynamics. The infinitesimal evolution operator is given 

as, 

t7(At) _ giAt(Lsystem + L NH C+LGLE) (27) 
rvj g*^ (LqlE+LnHc) gi&tLsystem g*"^ {^GLE+^N H c) 

(28) 

Where L is the Liouville operator. The integrators for MD can be obtained from the Trotter 
factorization of Liouville propagators. The corresponding Liouville operators for NH and 



GLE are given in detail in Refs. (|ll 2o|) respectively. The evolution of NH and GLE is 



updated at At/2 before and after the velocity- verlet routine for system evolution. As pointed 
out by Bussi and Parinellc 19 there is a significant drift in the conserved quantity of the GLE 
dynamics for strong friction coefficient (7 -1 ~ 10 fs). This is mainly due to the error 
introduced by the approximate integrator for GLE when jAt is not negligible. This error 
arises from the integration of the Hamilton equations and not from the friction itself. This 
problem can be solved by choosing a smaller time step for the simulation. In our simulation 
NGLE C (see Table. (JTJ) ) we use strong friction coefficient for zero point equilibration. For a 
stable conserved energy H' in Eq. ( TTTj) . we perform this simulation with a time step of 0.05 
fs. All the other simulations in this work are done with 0.5 fs time step. 

As described in section III Al we have two thermostats acting on the system and they 
compete to enforce their respective temperatures. The strength of the GLE thermostat 
varies with the frequency and is strongest for the modes in the range Uq ± Auo. The result 
of this competition is that the modes in the range ~ uq±Aoj are thermostated at an effective 
temperature T e // < Tnh- Depending on the strength of the friction coefficient (height 
of K(u) peak), the effective temperature of the specific modes can lie anywhere between 
< T(cj) e ff < T/vh- Temperature of all the other modes are equilibrated at T ~ T^h as the 
K(uj) ~ for u ^ (uq — Au, Uq + Au>). Hence we can control the effective temperature of a 
particular normal mode using NGLE dynamics. Our system consists of 256 water molecules 
with a density of 0.997 g/cm 3 modeled by the flexible force field model q-TIP4P/F 1 . This 
model has been used in many isotope effects studies both for water and ice^ 1 ^, and therefore 
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is an excellent model to evaluate the performance of our approach. We now demonstrate 
effects of damping of narrow range of modes using NGLE thermostats. This example will 
establish the spirit in which we intend to use NGLE thermostats. We use in house developed 
code for the force field and MD implementations. In Fig.Q]we show the vibrational density of 
states (obtained from the Fourier transform of the velocity autocorrelation function) for the 
q-TIP4P/F model. The three peaks corresponds to translation+rotations (400 — 1000 cm -1 ), 
bending (~ 1600 cm" 1 ) and stretching (~ 3600 cm -1 ) modes. We also plot the projected 
temperature of translation, rotation and intramolecular vibrational modes(see Fig. (ITTa))). 
Mode-projected temperatures is an important parameter to monitor NGLE action on the 
system. We calculate these projections by defining new molecular subspaces along the center 
of mass (for the translations), molecular main three moment of inertia axis (for the rotations) 
and the three vibrational normal modes of the isolated molecule. These projections act as a 
guide to tune NGLE parameters for designing a frequency dependent temperature control. 




5000 



at (cm l ) 



Figure 1: Spectral density plot of q-TIP4P/F liquid water equilibrated with NH thermostats at 300 
K 



For GLE implementation we base our code on codes developed by Ceriotti et al.—. We see 
that without the GLE action all the modes are equilibrated to the temperature set by NH 
thermostats. We code the memory kernels to overlap with the vibrational spectrum of the 
system and tailor them to our requirement. We select a memory profile of very narrow width 
sharply peaked at some frequency, uq which vanishes rapidly away from it. NH chains are set 
to 300 K. GLE is much more strongly coupled to the system than the NH chains at selective 
modes (oj ± Au). In the expression of memory kernel, Eq. (1201) . we set the peak position 
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Uq = 3600 cm -1 , peak width Au> = 5 cm -1 and the strength of coupling 7 _1 = 1000 ps (see 
Fig. [2]). We now extend this method to study how the structure of liquid water changes 
when modes are kept at different temperatures. 
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Figure 2: Spectral density plot of q-TIP4P/F liquid water at 300 K with NGLE, black line. Red 
line, delta-Like peak frequency depdendent friction profile. 



III. NGLE APPLIED TO Q-TIP4P/F WATER 

Using NGLE to simulate liquid water to non-equilibrium distribution of temperatures, we 
can analyze how the structure of q-T!P4P/F water depends on the temperature of individual 
modes. This way we can establish the influence of each individual mode in the structure 
of the water. Also one can study how the different dynamical modes are connected to each 
other. 

A. Weak damping of intramolecular modes 

We damp the intramolecular modes and study their influence on the overall structure of 
liquid water, using a damping profile as shown in Fig. ([3]). For this simulation (NGLE A) 
parameters are described in Tabled Projected temperatures plot (Fig [7(b)) show that the 
intramolecular modes are kept are relatively lower temperatures for the full simulation run. 
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Figure 3: Spectral density plot. The black line is for of q-TIP4P/F liquid water at 300 K, the 
solid red line is for the NGLE A (see Table [T]). The red dashed line shows the frequency dependent 
friction profile. 

We also plot the radial distribution function (rdf) for NGLE A. This rdf clearly shows a 
local softening of the first two 0-0 peaks. On the other hand the first peak of the O-H rdf 
as expected is much sharper. The higher order peaks, linked to the Hbond network have 
also become softer. A more frozen covalent bond results in a loss of structure of liquid water. 
This is a consequence of the previously mentioned anticorrelation. The Hbond is weakened 
by strengthening the O-H covalent bond. 



B. Weak damping of inter molecular modes 

In this section we analyze whether a more fluctuating O-H covalent bond induces a 
local structuring of liquid water. To answer this question we now damp the low energy 
intermolecular modes. We couple GLE to low energy modes with almost no coupling to 
intramolecular modes. The parameters for this simulation (NGLE B) are described in Tabled] 
The Results are shown in Figs. ( |5|7|) . 

The action of damping low energy modes results in relatively higher temperature of vibra- 
tional modes (see Fig. [7(c)). Consequently, we see more structure in the 0-0 rdf (see Fig. 
[6]). This establishes the well know anticorrelation 5 between inter and intramolecular modes 
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r(A) 

Figure 4: Radial distribution function plots of q-TIP4P/F liquid water for NVT simulation (black) 
at 300 K and NGLE damped OH stretching (red). Top: 0-0 rdf. Bottom: O-H rdfs, the inset 
shows a zoom into the first peak. 




1000 2000 3000 4000 5000 
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Figure 5: Spectral density plot. The black line is for of q-TIP4P/F liquid water at 300 K, the 
solid red line is for the NGLE B (see Table [I]). The red dashed line shows the frequency dependent 
friction profile. 

in liquid water. Hence we demonstrate that the strengthening of the water structure (by 
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more frozen Hbonds) implies a more delocalized 0-H intramolecular bond, i.e. a softening of 
the stretching OH vibration. In Table ([!]) we present normal mode-projected temperatures 
for translations, rotations and vibration modes. 
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Figure 6: Radial distribution function plots of q-TIP4P/F liquid water for NVT simulation (black) 
at 300 K and NGLE damped intermolecular stretching (red). Upper panel: 0-0 rdf. Lower panel 
: O-H rdfs, the inset shows a zoom into the first peak. 



C. Modes close to their zero point temperature 

None of the results shown in the previous sections were surprising, they are a confirmation 
of the anti-correlation effect. This effect is a manifestation of the strong anharmonicity of 
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Table I: Temperature distribution among individual modes 



T NH (K) 


Trans(K) 


Rot(K) 


Vib(K) 


(cm 1 ) 


Au (cm 1 ) 


7 1 (ps) 


Tqle (K) 
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Figure 7: Temperature of translational (black), rotational (red) and vibrational (green) modes in 
liquid water plotted as a function of time. This mode projected temperature is plotted for (a) 
NVT simulation at 300 K. (b) NGLE A simulation with damped intramolecular OH stretching . 
(c) NGLE B simulation with damped intermolecular modes, (d) NGLE C simulation with modes 
close to effective zero point temperature. See Table ([!]) for temperature values. 

the vibrational modes, which strongly couple to the rotational modes when Hbonds are 
formed. However, when all the normal modes are equilibrated at their corresponding zero 
point temperatures, the two opposite effects we just described should balance out. As shown 
by Habershon et a/.-, this competion results in an overall minimization of quantum effects on 
the structure of liquid water. To evaluate this using our method, we set the NH-temperature 
close to the zero point temperature of vibrational modes and damp intermolecular modes 
so that the effective temperature of the low energy modes are close to their effective zero 
point temperature. The shape of the tailored frequency dependent memory profile is shown 
in Fig. ([8]) . The advantage we have is that very few parameters are used to achieve this non 
equilibrium T distribution. 

We plot the rdfs of liquid water for this case (see Fig. ([9])) and compare them with 
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those obtained from a "classical" and a PIMD simulation, both of an identical system. By 
"classical" we mean a standard NVT MD simulation at T=300 K, with the use of a NH 
thermostat (without any GLE). The PIMD simulation was performed using the same method 
as described in Ref. (2j), using a polymer ring of P=32 beads. Fig. (jH]) clearly demonstrates 
that using NGLE (with NGLE-C parameters, see Table (JI|)) we can equilibrate modes to 
their average zero point temperature and the obtained rdfs are in good agreement with 
PIMD results. 

1. Vibrational spectrum 

Since NGLE is a deterministic thermostat, one of the advantages of this method is that 
dynamical properties are described with the same accuracy as in NH simulations. This al- 
lows us to study the changes on the vibrational spectrum induced by the new distribution 
of normal mode temperatures. In other words, it allows us to evaluate the intrinsic anahar- 
monicities of liquid water (which strongly depend on the underlying model) at the correct 
zero-point temperature distribution. 




Figure 8: Spectral density plots, decomposed into different vibrational contributions. Total spec- 
trum (top left panel). Projection onto translations (top right), rotations (bottom left) and vibra- 
tions (bottom right). Solid black line, NVT simulation at T=300 K. Red solid line, NGLE C (see 
Table (|T])) simulation. The red dashed line shows the frequency dependent friction profile used in 
NGLE C. 
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In order to accurately evaluate how each region of the vibrational spectrum changes, we 
have obtained the spectral density projected onto translational, rotational and vibrational 
modes separately. The procedure is straightforward once the cartesian coordinates are pro- 
jected onto the normal mode coordinates. The spectral density is obtained by computing 
the Fourier transform of the velocity-velocity autocorrelation function within each indepen- 
dent normal mode subspace. Results are presented in Fig. (jSJ). The figure shows both the 
changes to the total spectrum (top left panel), and to the translations (top right), rotations 
(bottom left) and vibrations (bottom right). 

The partition of the spectrum helps on identifying which are the modes that undergo 
major frequency shifts upon addition of quantum effects. The largest, and more interesting 
changes occur in the translational and rotational regions of the spectrum. While the lowest 
frequency translational motions remain unaffected, those modes with classical frequencies 
above ?^200 cm -1 are strongly blue shifted. These are modes associated to the Hbond 
stretching. The rotational peak also undergoes a large blue shift, although some of the 
modes also remain unaffected. The shifted rotational modes most likely correspond to 
those associated to Hbond bendings. In both cases, the shift is a consequence of the large 
temperature NGLE imposes on the intramolecular vibrational modes. Indeed, as seen in 
Table label "NGLE C" and in Fig. [7] (d), the net T of the vibrational modes is 2200 
K. This large blue shift confirms the strong coupling between inter and intra molecular 
modes, or a strong intrinsic anharmonicity. This shift also modifies their corresponding zero 
point temperature. To account for this zero point shift, our NGLE parameters are tuned to 
equilibrate intermolecular modes to 600 K (trans) 800 K (rotations) whereas intramolecular 
modes are equilibrated at 2200 K. 

The strong coupling between the different modes makes it difficult to fully determine the 
actual quantum vibrational spectrum of the model. How much a given model's vibrational 
spectrum should be modified upon consideration of NQE is an open question, not easy 
to address within a thermostated simulations^. Nonetheless, our scheme is particularly 
suited to extract dynamical information. First of all, the method ensure the absence of 
zero-point leakage, as seen in the conservation of mode-projected temperatures through the 
simulation length. Secondly, due to the suppression of GLE noise, we do not add any more 
diffusion to the system, other than that associated to the NH thermostat. A first principles 
MD simulation of liquid water and its vibrational spectrum under the action of our NGLE 
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thermostat will help to gain insight into the model-dependence of the NGLE spectrum. 



IV. ANALYSIS AND DISCUSSION 

In this section we present analysis of our results of selective mode thermostating of liquid 
water. Based on our results we are in a position to provide deeper insight into the role of 
individual modes contributing to the overall structure of liquid water. 

A. Emergence of water structure from stretching modes 

The results for zero point simulation using NGLE can also be interpreted from a different 
view point. In this section we compare classical simulation done at 600 K to our effective zero 
point simulation done at 600 K (translations), 800 K (rotations) and 2200 K (vibrations) 
("NGLE C" in Table ©). 



Fig. (ITUj) shows that for a classical simulation equilibrated at 600 K the liquid water structure 
is completely lost, as expected for such a high temperature. Indeed at this T the liquid is at 
a supercritical state, because we do not allow the volume of our simulation cell to change. 
The 0-0 rdf in red is a liquid water structure that emerges completely due to the higher 
temperature imposed on the stretching modes. This comparison clearly points out that zero 
point temperature of O-H stretching modes plays a dominant role in the structuring of liquid 
water. 

B. Role of individual modes from PIMD simulations. 

So far in this work we have demonstrated how quantum effects corresponding to individual 
modes effect the overall structure of liquid water. In this section we try to extract this 
information from PIMD simulations. Recent work of Habershon et a/.- established the 
idea of competing quantum effects by understanding the role of intramolecular stretching 
modes. They observe reduced quantum to classical ratio of diffusion coefficient when the 
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Figure 9: Radial distribution plots. Top: 0-0 rdf. Middle O-H rdf. Bottom H-H rdf. Black line: 
NVT simulation at 300 K. Red line: Zero-point NGLE simulation (NGLE C, in Table ©). Blue 
line: PIMD simulation with 32 beads. 

O-H stretching is allowed. This is due to the anharmonicity of the stretching mode, which 
couples it to the rotational and translational modes. This is in agreement with our results. 
However in this work we aim to understand this effect in terms of the zero point temperature 
of competing modes. In PIMD simulations, we map the zero point temperature on individual 
modes to the number of beads. The quantum limit is achieved with P — > oo, P being the 
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Figure 10: 0-0 Radial distribution function. Red line: NGLE C (see Table (J])) with a temperature 
distribution of: 600 K translations, 800 K rotations and 2200 K vibrations. Black line: NVT 
simulation equilibrated at 600 K. 

number of beads employed for the discretization of the path integral. A PIMD simulation 
with finite P at temperature T implies a high-temperature approximation, i.e. the partition 
function of the system is considered to be classical at a higher temperature given by the 
product PT. Usually the value of P is chosen so that the product PT is several times larger 
than the zero point temperature of the highest frequency. 

We study the structure of liquid water as a function of P. We plot the rdf obtained from 
PIMD for 0-0 and O-H pairs, for a 256 molecules water simulation of 250 ps length and 
density 0.997 g/cm 3 (see Fig. [TT]) . 



Based on the plots in Fig. ( [111 we can analyze how the zero point effects of each mode 
influences the structure of liquid water. To understand the structure we introduce an order 
parameter which is the ratio of the first maxima to the first minima of the 0-0 rdf. The 
increasing value of this ratio is related to the increasing structure of liquid water. We plot 
this ratio as a function of P. For the O-H rdf we plot the height of the first peak as a 
function of P (see Fig. (1121)). 

In Fig. ( I12j) we observe that for P < 10, i.e. for PT < 3000K goo(r) looses structure 
compared to the classical (P=l bead) case. This can be attributed to the fact that using 
only a few beads, the zero point temperature of intermolecular modes is well captured and 
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Figure 11: Radial distribution function plots for PIMD simulations as a function of the number of 
beads. Top: 0-0 rdf. The two insets are zoomed into the first peak and the first minima. Bottom: 
O-H rdf 

they tend to unstructure the liquid. However for small P the PIMD simulation is not able to 
describe the zero point vibration of the O-H stretching modes. We have already mentioned 
that the temperature PT must be several times larger that the zero point temperature for 
a reasonable approximation of the quantum limit in a PIMD simulation.. As we further 
increase no. of beads we see that goo(r) starts gaining structure. This gain in structure is 
due to the better and better description of zero point effects of intramolecular modes. This 
gain in structure saturates when the intramolecular modes are kept close to their zero point 
temperature. Hence using P as a control parameter for quantum effects, we can understand 
the competition between intra and intermolecular modes. We also see a crossover from 
unstructuring to structuring in goo(r) (Fig. (TTTT) ) as we increase the no. of beads. A model 
with too little anticorrelation as the one used here£, will not structure above the classical 
level, independently of the number of beads used. 
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Figure 12: Order parameter as a function of the number of beads in PIMD simulations. Top: ratio 
of the height of first maxima to first minima. Bottom: height of first peak of the O-H rdf. The 
points are the actual data and the solid line is a fourth order interpolation. 

V. CONCLUSIONS 

In this work, we have combined NH and GLE thermostats to create a powerful selective 
mode thermostating scheme, that we coin NGLE. Using NGLE we were able to equilibrate 
vibrational modes to different temperatures, ensuring the maintenance of the NGLE-imposed 
non-equilibrium temperature distribution for any simulation time length. The GLE noise is 
suppressed by setting the GLE temperature close to zero, making the thermostat dynamics 
deterministic. We have applied NGLE to a flexible force field model of water (q-TIP4P/F). 
Our results show that the structure of liquid water changes when intramolecular modes are 
equilibrated to a different temperature than that of intermolecular modes. We showed that 
equilibrating vibrations at slightly lower temperature than other modes results in the loss of 
structure in liquid water. Conversely water is more structured if we equilibrate vibrations 
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to slightly higher temperature. This simple exercise verifies the well known anticorrelation 
effect 5 in terms of the temperature of individual modes. Simply changing NGLE parame- 
ters we can set all modes to their corresponding zero point temperature. Equilibrating the 
intramolecular modes close to their zero point temperature (~ 2200 K) results in blue shift 
in the peak of the intermolecular (translation+rotation) modes. This indicates the large 
intrinsic anharmonicity of the vibrational modes in liquid water. Our zero point estimates 
of translations and rotations prove to be correct as NGLE simulations reproduce O-O, O-H 
and H-H rdf obtained from PIMD simulation with 32 beads. Finally, we have analyzed 
selective mode behavior as a function of their zero point temperature in PIMD. To this aim 
we showed that by a systematic increase of the number of beads P in PIMD simulations, the 
different vibrational modes of water can be successively tuned from a classical to a quantum 
limit. The height of first peak of the O-H rdf continuously decreased as a function of P. 
This effect is similar to adding more zero point temperature on vibrations. For P < 10, 
there is drastic softening of long range structure as seen in the 0-0 rdf. This is the case 
when the temperature of low energy modes are close to their zero point temperature but 
the vibration temperature is still far away from its zero point. For P > 10 we observe an 
increase of structure in the 0-0 rdf that asymptotically saturates with P. This is a conse- 
quence of intramolecular modes equilibrated to their zero point temperature. In summary, 
we successfully introduced a selective mode thermostating scheme (NGLE), which can be 
easily tuned to include quantum zero point effects, producing results in good agreement 
with PIMD. We used it to study how the structure of liquid water responds to different 
temperatures, unveiling the existence of large intrinsic anharmonicities in the vibrational 
modes of liquid water. 
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